#################################
### Revisiting White Backlash ###
###     Replication File      ###
###        TESS Data           ###
#################################
library(foreign)
library(car)
library(MASS)
library(sandwich)
library(AER)
library(texreg)
library(multcomp)
library(stargazer)
library(pwr)

# Load data
Data <- read.spss("Tess203_Piston_Client.sav", 
                  use.value.labels=TRUE, 
                  to.data.frame = TRUE)

##################
### Clean data ###
##################

# Racial resentment:
# Irish, Itialian, Jewish and many other minorities overcame prejudice and worked
# their way up. Blacks should do the same without any special favors.
Data$Q9_1 <- recode(Data$Q9_1,
                    "'Strongly Agree' = 2;
                    'Agree' = 1;
                    'Neither Agree nor Disagree' = 0;
                    'Disagree' = -1;
                    'Strongly Disagree' = -2;
                    'Refused' = NA",
                    as.factor.result = F)
# Generations of slavery and discrimination have created conditions that make it
# difficult for blacks to work their way out of the lower class
Data$Q9_2 <- recode(Data$Q9_2,
                    "'Strongly Agree' = -2;
                    'Agree' = -1;
                    'Neither Agree nor Disagree' = 0;
                    'Disagree' = 1;
                    'Strongly Disagree' = 2;
                    'Refused' = NA",
                    as.factor.result = F)
# Over the past few years, blacks have gotten less than they deserve
Data$Q9_3 <- recode(Data$Q9_3,
                    "'Strongly Agree' = -2;
                    'Agree' = -1;
                    'Neither Agree nor Disagree' = 0;
                    'Disagree' = 1;
                    'Strongly Disagree' = 2;
                    'Refused' = NA",
                    as.factor.result = F)
# It's really a matter of some people not trying hard enough; if blacks would 
# only try harder they would be just as well off as whites
Data$Q9_4 <- recode(Data$Q9_4,
                    "'Strongly Agree' = 2;
                    'Agree' = 1;
                    'Neither Agree nor Disagree' = 0;
                    'Disagree' = -1;
                    'Strongly Disagree' = -2;
                    'Refused' = NA",
                    as.factor.result = F)
# Racial resentment: low -> high; -8 -> 8
Data$Resentment <- Data$Q9_1 + Data$Q9_2 + Data$Q9_3 + Data$Q9_4
# PID
Data$party7 <- recode(Data$party7,
                      "'Strong Democrat' = 7;
                      'Not Strong Democrat' = 6;
                      'Leans Democrat' = 5;
                      'Undecided/Independent/Other' = 4;
                      'Leans Republican' = 3;
                      'Not Strong Republican' = 2;
                      'Strong Republican' = 1;
                      'Refused' = NA",
                      as.factor.result = F)
# Ideology
Data$ideo <- recode(Data$ideo,
                    "'Extremely liberal' = 7;
                    'Liberal' = 6;
                    'Slightly liberal' = 5;
                    'Moderate, middle of the road' = 4;
                    'Slightly conservative' = 3;
                    'Conservative' = 2;
                    'Extremely conservative' = 1;
                    'Refused' = NA",
                    as.factor.result = F)
# Education
Data$PPEDUC <- recode(Data$PPEDUC,
                      "'Not asked' = NA;
                      'REFUSED' = NA;
                      'No formal education' = 1;
                      '1st, 2nd, 3rd, or 4th grade' = 1;
                      '5th or 6th grade' = 1;
                      '7th or 8th grade' = 1;
                      '9th grade' = 1;
                      '10th grade' = 1;
                      '11th grade' = 1;
                      '12th grade NO DIPLOMA' = 1;
                      'HIGH SCHOOL GRADUATE - high school DIPLOMA or the equivalent (GED)' = 2;
                      'Some college, no degree' = 3;
                      'Associate degree' = 4;
                      'Bachelors degree' = 5;
                      'Masters degree' = 6;
                      'Professional or Doctorate degree' = 7",
                      as.factor.result = F)
# Gender
Data$PPGENDER <- recode(Data$PPGENDER,
                        "'Not asked' = NA;
                        'REFUSED' = NA;
                        'Male' = 1;
                        'Female' = 0",
                        as.factor.result = F)
# Income
Data$PPINCIMP <- recode(Data$PPINCIMP,
                        "'Not asked' = NA;
                        'REFUSED' = NA;
                        'Less than $5,000' = 1;
                        '$5,000 to $7,499' = 1;
                        '$7,500 to $9,999' = 1;
                        '7th or 8th grade' = 1;
                        '$10,000 to $12,499' = 1;  
                        '$12,500 to $14,999' = 1;
                        '$15,000 to $19,999' = 1; 
                        '$20,000 to $24,999' = 2;
                        '$25,000 to $29,999' = 2;  
                        '$30,000 to $34,999' = 2;   
                        '$35,000 to $39,999' = 3;  
                        '$40,000 to $49,999' = 3;   
                        '$50,000 to $59,999' = 4;
                        '$60,000 to $74,999' = 4;   
                        '$75,000 to $84,999' = 5;
                        '$85,000 to $99,999' = 5;  
                        '$100,000 to $124,999' = 6;
                        '$125,000 to $149,999' = 6; 
                        '$150,000 to $174,999' = 6;
                        '$175,000 or more' = 6",
                        as.factor.result = F)
# Age
Data$PPAGE <- ifelse(Data$PPAGE > 79, 7, 
                     ifelse(Data$PPAGE > 69, 6,
                            ifelse(Data$PPAGE > 59, 5,
                                   ifelse(Data$PPAGE > 49, 4,
                                          ifelse(Data$PPAGE > 39, 3,
                                                 ifelse(Data$PPAGE > 29, 2, 1))))))
# Received treatment frame
Data$Treatment <- is.na(Data$Q8A)

# Support death penalty
Data$SupportDP <- ifelse(is.na(Data$Q8A), as.character(Data$Q8B), as.character(Data$Q8A))
Data$SupportDP <- recode(Data$SupportDP,
                         "'Strongly favor' = 7; 
                         'Somewhat favor' = 6;
                         'Slightly favor' = 5; 
                         'Neither favor nor oppose' = 4;
                         'Slightly oppose' = 3;
                         'Somewhat oppose' = 2;
                         'Strongly oppose' = 1;
                         'Refused' = NA",
                         as.factor.result = F)

# Binary favor dp
Data$Favor <- ifelse(Data$SupportDP < 4, 0, ifelse(Data$SupportDP == 4, NA, 1))

# Remove other treated units
Control <- Data[Data$XTESS203 == 'control', c('SupportDP', 'Treatment', 'Favor',
                                              'Resentment', 'party7', 'ideo', 
                                              'PPEDUC', 'PPETHM', 'PPGENDER',
                                              'PPINCIMP', 'PPAGE')]

# Subset whites (an unnecessary step)
Whites <- Control[Control$PPETHM == 'White, Non-Hispanic', ]


##############################
### Descriptive Statistics ###
###         pg. 11          ###
##############################

### --- Gender --- ###
length(which(Whites$PPGENDER == 0))/length(Whites$PPGENDER) # female
#[1] 0.4977307

### --- Age --- ###
# Age 1=18-29, 2=30-39, 3=40-49, 4=50-59, 5=60-69, 6=70-79, 7=80+
summary(as.factor(Whites$PPAGE))
for(i in 1:7){
  print(length(which(as.factor(Whites$PPAGE) == i))/length(as.factor(Whites$PPAGE)))
}

# [1] 0.1270802
# [1] 0.1134644 +
# [1] 0.1210287 +
# [1] 0.2329803 = 0.4674734
# [1] 0.2481089 + 
# [1] 0.1255673 +
# [1] 0.03177005 = .4054463

### --- Education --- ###
# Education: 1=No HS, 2=HS, 3=Some college, 4=Associate degree, 5=Bachelors degree, 6=Masters degree, 7=Doctoral degree

for(i in 1:7){
  print(length(which(Whites$PPEDUC== i))/length(Whites$PPEDUC))
}

# [1] 0.05143722 +
# [1] 0.311649 +
# [1] 0.1754917 = .5385779
# [1] 0.09984871 +
# [1] 0.2133132 = .3131619
# [1] 0.1134644 +
# [1] 0.03479576 = .1482602

### --- Income --- ###
#Income 1=<20K, 2=20-35, 3=35-50, 4=51-75, 5=75-100, 6=100+, 7=NA
for(i in 1:6){
  print(length(which(Whites$PPINCIMP == i))/length(Whites$PPINCIMP))
}

# [1] 0.1164902
# [1] 0.1225416
# [1] 0.1270802
# [1] 0.1633888
# [1] 0.1422088
# [1] 0.3282905

### --- Ideology --- ###
# Ideology: Conservative -> Liberal
for(i in 1:7){
  print(length(which(Whites$ideo == i))/length(Whites$ideo))
}

# [1] 0.04992436 - strong cons.
# [1] 0.2239032
# [1] 0.1149773 cumulatively = .38888049
# [1] 0.3358548
# [1] 0.1059002 cumulatively = .2647505
# [1] 0.1270802
# [1] 0.03177005 - strong lib.

### --- PID --- ###
# Rep -> Dem
# PID 7: Others coded in the Independent/Category
for(i in 1:7){
  print(length(which(Whites$party7 == i))/length(Whites$party7))
}

# [1] 0.1664145 - strong R
# [1] 0.128593
# [1] 0.2133132 cumulatively = 0.5083207
# [1] 0.03479576
# [1] 0.1815431 cumulatively = 0.4568835
# [1] 0.1164902
# [1] 0.1588502 - strong D


###################################
### Ordered Probit: TESS Result ###
###        Table 4, p. 13        ###
###################################

OP <- polr(as.factor(SupportDP) ~ Treatment, data = Whites,
           method = 'probit', Hess = T)

# result

coeftest(OP)

texreg(OP, stars = numeric(0), include.thresholds = T, digits = 3, 
       custom.coef.names = c('Race frame',
                             'Strongly Oppose $|$ Somewhat Oppose',
                             'Somewhat Oppose $|$ Slightly Oppose',
                             'Slightly Oppose $|$ Neither Favor nor Oppose',
                             'Neither Favor nor Oppose $|$ Slightly Favor',
                             'Slightly Favor $|$ Somewhat Favor',
                             'Somewhat Favor $|$ Strongly Favor'))


##############################
### Bar Chart of % Support ###
###     Figure 2, p.14     ###
##############################

### --- Barchart of Percentage Differences --- ###
# t-test
ttest <- t.test(x = Whites[Whites$Treatment == FALSE, 'Favor'], 
                y = Whites[Whites$Treatment == TRUE, 'Favor'])
ttest$estimate[2] - ttest$estimate[1]

# vector of point estimates
PointEstimates <- c(ttest$estimate)*100

# confidence intervals for each treatment
Conf1 <-  t.test(Whites[Whites$Treatment == FALSE, 'Favor'])$conf.int
Conf2 <-  t.test(Whites[Whites$Treatment == TRUE, 'Favor'])$conf.int

# vector of confidence intervals
ConfInts <- list(Conf1, Conf2)

# create barplot
barplot(PointEstimates, ylim = c(0, 100), space = c(0),
        xlab = '', ylab = 'Percent favoring death penalty', axes = F, 
        names.arg = F, col = c('white', 'grey65'), 
        main = '')
# add y-axis
axis(2, las = 1)
# add legend
legend('top', bty = 'n', title = 'Frame', cex = 0.65,
       legend = c('None', 'Race'), fill = c('white', 'grey65'))
# add confidence bars
for (i in 1:2) {
  segments(x0 = c(0, 1)[i]+0.5, x1 = c(0, 1)[i]+0.5, 
           y0 = ConfInts[[i]][1]*100, y1 = ConfInts[[i]][2]*100)
  segments(x0 = c(0, 1)[i]+0.45, x1 = c(0, 1)[i]+0.55,
           y0 = ConfInts[[i]][1]*100, y1 = ConfInts[[i]][1]*100)
  segments(x0 = c(0, 1)[i]+0.45, x1 = c(0, 1)[i]+0.55,
           y0 = ConfInts[[i]][2]*100, y1 = ConfInts[[i]][2]*100)
}
# add text
text(0.5,10, 'N=')
text(0.5,5, length(Whites$Treatment)-sum(Whites$Treatment))
text(1.5,5, sum(Whites$Treatment))


####################
### OLS: TESS    ###
###   Table S2   ###
####################

OLS <- lm(SupportDP ~ Treatment, data = Whites)

# results
summary(OLS)

# robust SEs
sqrt(diag(vcovHC(OLS)))


###########################################
### Heterogeneous Treatment: TESS Data  ###
###                Table S12             ###
###########################################

# Republican regression
Whites$GOP <- Whites$party7 < 4
GOPOP <- polr(as.factor(SupportDP) ~ GOP*Treatment,
              data = Whites, 
              method = 'probit')
summary(GOPOP)
coeftest(GOPOP)
MEGOP <- glht(GOPOP, linfct = c("TreatmentTRUE + GOPTRUE:TreatmentTRUE = 0"))
summary(MEGOP)
# power = 0.543
pwr.t2n.test(n1 = nrow(Whites[Whites$GOP == T & Whites$Treatment == T, ]),
             n2 = nrow(Whites[Whites$GOP == T & Whites$Treatment == F, ]),
             d = 0.02718/0.12012,
             sig.level = 0.05,
             alternative = 'two.sided')

# Conservative regression
Whites$Con <- Whites$ideo < 4
ConOP <- polr(as.factor(SupportDP) ~ Con*Treatment,
              data = Whites, 
              method = 'probit')
summary(ConOP)
coeftest(ConOP)
MECon <- glht(ConOP, linfct = c("TreatmentTRUE + ConTRUE:TreatmentTRUE = 0"))
summary(MECon)
# power = 1
pwr.t2n.test(n1 = nrow(Whites[Whites$Con == T & Whites$Treatment == T, ]),
             n2 = nrow(Whites[Whites$Con == T & Whites$Treatment == F, ]),
             d = 0.1327/0.1372,
             sig.level = 0.05,
             alternative = 'two.sided')

# Low education regression
Whites$LowEdu <- Whites$PPEDUC < 4
LowEduOP <- polr(as.factor(SupportDP) ~ LowEdu*Treatment,
                 data = Whites, 
                 method = 'probit')
summary(LowEduOP)
coeftest(LowEduOP)
MELowEdu <- glht(LowEduOP, linfct = c("TreatmentTRUE + LowEduTRUE:TreatmentTRUE = 0"))
summary(MELowEdu)
# power = 1
pwr.t2n.test(n1 = nrow(Whites[Whites$LowEdu == T & Whites$Treatment == T, ]),
             n2 = nrow(Whites[Whites$LowEdu == T & Whites$Treatment == F, ]),
             d = 0.2266/0.1163,
             sig.level = 0.05,
             alternative = 'two.sided')

# High resentment regression
Whites$HiRes <- Whites$Resentment > 0
HiResOP <- polr(as.factor(SupportDP) ~ HiRes*Treatment,
                data = Whites, 
                method = 'probit')
summary(HiResOP)
coeftest(HiResOP)
MEHiRes <- glht(HiResOP, linfct = c("TreatmentTRUE + HiResTRUE:TreatmentTRUE = 0"))
summary(MEHiRes)
# power = 0.0729
pwr.t2n.test(n1 = nrow(Whites[Whites$HiRes == T & Whites$Treatment == T, ]),
             n2 = nrow(Whites[Whites$HiRes == T & Whites$Treatment == F, ]),
             d = 0.004659/0.107956,
             sig.level = 0.05,
             alternative = 'two.sided')

# Age regression
Whites$Young <- ifelse(Whites$PPAGE == 1, 1, 0)
Whites$Middle <- ifelse(Whites$PPAGE > 1 & Whites$PPAGE < 5, 1, 0)
Whites$Old <- ifelse(Whites$PPAGE > 4, 1, 0)
YoungOP <- polr(as.factor(SupportDP) ~ Young*Treatment,
                data = Whites, 
                method = 'probit')
MiddleOP <- polr(as.factor(SupportDP) ~ Middle*Treatment,
                 data = Whites, 
                 method = 'probit')
OldOP <- polr(as.factor(SupportDP) ~ Old*Treatment,
              data = Whites, 
              method = 'probit')
summary(YoungOP)
coeftest(YoungOP)
MEAge1 <- glht(YoungOP, linfct = c("TreatmentTRUE + Young:TreatmentTRUE = 0"))
summary(MEAge1)
# power = 0.588
pwr.t2n.test(n1 = nrow(Whites[Whites$Young == 1 & Whites$Treatment == T, ]),
             n2 = nrow(Whites[Whites$Young == 1 & Whites$Treatment == F, ]),
             d = 0.1115/0.2308,
             sig.level = 0.05,
             alternative = 'two.sided')

MEAge2 <- glht(MiddleOP, linfct = c("TreatmentTRUE + Middle:TreatmentTRUE = 0"))
summary(MEAge2)
# power = 0.946
pwr.t2n.test(n1 = nrow(Whites[Whites$Middle == 1 & Whites$Treatment == T, ]),
             n2 = nrow(Whites[Whites$Middle == 1 & Whites$Treatment == F, ]),
             d = 0.04969/0.12193,
             sig.level = 0.05,
             alternative = 'two.sided')

MEAge3 <- glht(OldOP, linfct = c("TreatmentTRUE + Old:TreatmentTRUE = 0"))
summary(MEAge3)
# power = 1
pwr.t2n.test(n1 = nrow(Whites[Whites$Old == 1 & Whites$Treatment == T, ]),
             n2 = nrow(Whites[Whites$Old == 1 & Whites$Treatment == F, ]),
             d = 0.2025/0.1318,
             sig.level = 0.05,
             alternative = 'two.sided')

texreg(list(GOPOP, ConOP, LowEduOP, HiResOP, YoungOP, MiddleOP, OldOP), stars = numeric(0), include.thresholds = T, digits = 3, 
       custom.coef.names = c('Republican',
                             'Race frame',
                             'Republican $\\times$ race frame',
                             'Strongly Oppose $|$ Somewhat Oppose',
                             'Somewhat Oppose $|$ Slightly Oppose',
                             'Slightly Oppose $|$ Neither Favor nor Oppose',
                             'Neither Favor nor Oppose $|$ Slightly Favor',
                             'Slightly Favor $|$ Somewhat Favor',
                             'Somewhat Favor $|$ Strongly Favor',
                             'Conservative', 
                             'Conservative $\\times$ race frame',
                             'Low education', 
                             'Low education $\\times$ race frame',
                             'High racial resentment', 
                             'High racial resentment $\\times$ race frame',
                             'Age 18-29',
                             'Age 18-29 $\\times$ race frame',
                             'Age 30-59',
                             'Age 30-59 $\\times$ race frame',
                             'Age 60+',
                             'Age 60+ $\\times$ race frame'))


######################
### Balance Checks ###
###   Table S14    ###
######################

PIDBal <- lm(party7 ~ Treatment, data = Whites)
IdeolBal <- lm(ideo ~ Treatment, data = Whites)
EducBal <- lm(PPEDUC ~ Treatment, data = Whites)
GenderBal <- lm(PPGENDER ~ Treatment, data = Whites)
IncomeBal <- lm(PPINCIMP ~ Treatment, data = Whites)
AgeBal <- lm(PPAGE ~ Treatment, data = Whites)

summary(PIDBal)
summary(IdeolBal)
summary(EducBal)
summary(GenderBal)
summary(IncomeBal)
summary(AgeBal)


####################################
### Bar Chart of Ordinal Support ###
###           Figure S2          ###
####################################

# t-test
ttestOrd <- t.test(x = Whites[Whites$Treatment == FALSE, 'SupportDP'], 
                   y = Whites[Whites$Treatment == TRUE, 'SupportDP'])
ttestOrd$estimate[2] - ttestOrd$estimate[1]

# vector of point estimates
PointEstimatesOrd <- c(ttestOrd$estimate)

# confidence intervals for each treatment
OrdConf1 <-  t.test(Whites[Whites$Treatment == FALSE, 'SupportDP'])$conf.int
OrdConf2 <-  t.test(Whites[Whites$Treatment == TRUE, 'SupportDP'])$conf.int

# vector of confidence intervals
ConfIntsOrd <- list(OrdConf1, OrdConf2)

# create barplot
barplot(PointEstimatesOrd, ylim = c(1, 7), space = c(0), xpd = F,
        xlab = '', ylab = 'Average death penalty support', axes = F, names.arg = F, 
        col = c('white', 'grey65'), 
        main = '')
# add y-axis
axis(2, las = 1)
# add legend
legend('top', bty = 'n', title = 'Frame', cex = 0.65,
       legend = c('None', 'Race'), fill = c('white', 'grey65'))
# add confidence bars
for (i in 1:2) {
  segments(x0 = c(0, 1)[i]+0.5, x1 = c(0, 1)[i]+0.5, 
           y0 = ConfIntsOrd[[i]][1], y1 = ConfIntsOrd[[i]][2])
  segments(x0 = c(0, 1)[i]+0.45, x1 = c(0, 1)[i]+0.55,
           y0 = ConfIntsOrd[[i]][1], y1 = ConfIntsOrd[[i]][1])
  segments(x0 = c(0, 1)[i]+0.45, x1 = c(0, 1)[i]+0.55,
           y0 = ConfIntsOrd[[i]][2], y1 = ConfIntsOrd[[i]][2])
}
segments(x0=0, x1=2, y0=1, y1=1, lwd=2)
# add text
text(0.5,1.45, 'N=')
text(0.5,1.15, length(Whites$Treatment)-sum(Whites$Treatment))
text(1.5,1.15, sum(Whites$Treatment))
